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Abstract 

We give a quick and dirty, but reasonably safe, algorithm for the minimization of a 
convex quadratic function under convex quadratic constraints. The algorithm minimizes 
the Lagrangian dual by using a safeguarded Newton method with non-negativity 
constraints. 
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1 Introduction 

In this note we give an algorithm for minimizing a quadratic / 0 over all x for which f s (x) < 0, 
where the f s with s = 1, • • • ,p are also quadratics. We assume the f s are convex quadratics 
for all s — 0, • • • ,p. Thus we have f s (x) = c s + b' s x + | x'A s x, with A s positive semi-definite 
for all s. In addition we assume, without loss of generality, that Aq is non-singular (and thus 
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positive definite). Finally, we assuc the Slater condition is satisfied, i.e. there is an x such 
that f s (x) < 0 for s — 1, • • • ,p. 

Quadratic programming with quadratic constraints (QPQC) has been studied in great detail, 
both for the convex and the much more complicated non-convex case. The preferred algorithm 
translates the problem into a second-order cone programming (SOCP) problem, which is then 
solved by interior point methods. This is the method used in the R packages DWD (Huang et 
al. (2011)), CLSOCP (Rudy (2011)), and cccp (Pfaff (2015)). 

We go a somewhat different route, more direct, but undoubtedly less efficient. In the appendix 
we collect some general expressions for the first and second derivatives of a marginal function 
and apply them to the Lagrangian of a constrained optimization problem. In the body of 
the paper we apply these formulas to the Langrangian dual of the QPQC problem, using a 
straightforward version of the Newton method with non-negativity constraints. 


2 QPQC 


Suppose the f s are convex quadratics for all s = (),••• ,p. We have f s {x) = c s + b' s x + \x'A s x, 
with A s positive semi-definite for all s. In addition we assume, without loss of generality, 
that A 0 is non-singular (and thus positive definite). It follows that the Lagrangian is 


p p i r p 

g (x, y ) = (c 0 + y»°s) + ( b o + ysb s )'x + -zx' < a + ys A s 

S=1 S=1 Z { .5=1 


x, 


and thus 


( p \ - 1 p 

x(y) = ~ < A) + Vs A s \ {bo + ysbs), ( 1 ) 

l S=1 J 5=1 

and 


p p f p ) _1 p 

h{y ) = (co + ys°s) - + ysb a y Mo + y* A * \ ( b o + J2 ysbs)- (2) 

S=1 Z 5 = 1 l 5=1 J 5=1 

We could use the explicit formula (2) to compute derivatives of h, but instead we apply the 
results of the appendix. 

In the QPQC case we have, from (9), 


Vh(y) = F{x{y)) 


For the second derivatives we need 


M x (y)) 


p 

g(x, y) = A 0 + J2 ys A s, 

S= 1 
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and 


Thus, from (10), 


VF{x ) 


b[ + x'Ai 
b' v + x'Ap 


{V 2 Ky)} st = 


~(A s x(y) + b s y 


Ao + ^2 y s A s 


s= 1 


(A t x(y) + b t ) 


We now know how to compute first and second derivatives of the dual function. In step k of 
the iterative process we use the quadratic Taylor approximation at the current value y^ of y. 

h(y) = % ( ‘>) + (y- + \(y - y {t) )'V\y^)(y - 

We then maximize this over y > 0 the find p, using the pnnqp function from the lsei 
package (Wang, Lawson, and Hanson (2015)). 

Finally we stabilize Newton’s method by computing 

r/( fc+1 ) = argmin h(y^ + A (y^ — y^)). 

0<A<1 

For this step-size or relaxation computation we use optimize from base R. But first we check 
if (— y^)'Vh(y^) > 0, in which case we choose step size equal to one. 

3 Example 

We give a small and artificial example. To make sure the Slater condition is satisfied we 
choose c s < 0 for s — 1, • • • ,p, which guarantees that f s ( 0) < 0. The example has quadratics 
of three variables and five quadratic constraints. 

set.seed(54321) 
a <- array (0, c(3, 3, 6)) 
for (j in 1:6) { 

a[, , j] <- crossprod (matrix (rnorm (300), 100, 3)) / 100 

} 

b <- matrix (rnorm (18), 3, 6) 
c <- -abs (rnorm (6)) 

print (qpqc (c(0, 0, 0, 0, 0), a, b, c, eps = le-10, verbose = TRUE)) 


## 

Iteration: 

1 

fold 

-Inf 

fnew: 

-2.69125228 

step: 

1.00000000 

## 

Iteration: 

2 

fold 

-2.69125228 

fnew: 

-2.65100837 

step: 

1.00000000 

## 

Iteration: 

3 

fold 

-2.65100837 

fnew: 

-2.65032456 

step: 

1.00000000 

## 

Iteration: 

4 

fold 

-2.65032456 

fnew: 

-2.65032433 

step: 

1.00000000 

## 

Iteration: 

5 

fold 

-2.65032433 

fnew: 

-2.65032433 

step: 

1.00000000 
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## $h 

## [1] -2.650324329 
## 

## $f 

## [1] -2.650324329 
## 

## $xmin 

## [1] 0.6424151506 -1.6326182487 0.2190791799 

## 

## $multipliers 

## [1] 0.0000000000 0.0000000000 0.0000000000 0.1040112458 0.0000000000 
## 

## Sconstraints 

## [1] -7.755080598e-01 -4.529990503e+00 -9.397132870e-01 1.047517628e-ll 

## [5] -2.269181794e+00 

4 Appendix: Derivatives 

In this appendix we collect some formulas for first and second derivatives of minima and 
minimizers in constrained problems. It is just convenient to have them in a single place, and 
we can apply them in the body of the paper. 

4.1 Marginal Function 

Suppose g : M n (8) —» M, and define 


%) := min g(x,y) = g(x(y),y), 


( 3 ) 


with 


x(y) := axgmin g(x,y) 7 (4) 

X 

which we assume to be unique for each y. 

Now assume two times continuous differentiability of g. Then x(y) satisfies 

g{x{y),y) = 0, 

and thus 

'^ug(x(y),y)Vx(y)+V 12 g(x(y),y) = 0. 

It follows that 

Vx(y) = -V^g(x(y),y)V 12 g(x(y),y). (5) 
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Now 


Vh(y) = V 1 g(x(y),y)Vx(y) + V 2 g(x(y),y) = V 2 g(x(y),y), (6) 


and 


V 2 h(y) = V 12 g(x(y),y)Vx(y) +V 22 g(x(y),y) = 

= V 22 g(x(y),y) - V 21 g(x(y),y)V^g(x(y),y)V l 2 g(x(y),y). (7) 

4.2 Lagrangians 

We specialize the results of the previous section to Langrangian duals. For the problem of 
minimizing / : M n —> M over x such that f s (x ) < 0 for s — 1, • • • ,p , where f s : M n — > M. We 
can also write the inequality constraints compactly as F(x) < 0, where F : M n (8) — > M. 
The Lagrangian is 


g(x,y) = fo(x) + y'F(x). 


( 8 ) 


The primal problem is 


min max q(x,y) = min 

x ^>0 x 


The dual problem is 

max ruin g(x 1 y) = nra xh(y). 

y >0 x y> 0 

For any feasible x and y we have weak duality h(y) < fo(x). If the f s for s — 0, • • • ,p are 
convex and the Slater condition is satisfied we have equal optimal values for the primal and 
dual problems. Moreover if y solves the dual problen then x(y) solves the primal problem. 

Thus, from (6), 


f 0 (x) if F(x) < 0, 
+oo otherwise 


= min{/ 0 (x) | F(x) < 0}. 


Vh(y) = V 2 g(x(y),y) = F(x(y)), (9) 

and, from (7), 

V 2 h(y) = -V 21 g(x(y),y)V^ 1 1 g(x(y),y)V 12 g(x(y),y). 

Also 

£>11 g(x(y),y) = V 2 f 0 (x(y)) + Vs^ 2 f s (x(y)), 

s= 1 

and 

V 21 g(x(y),y) = VF(x(y)\ 

so that 
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-1 


( 10 ) 


V 2 h(v) = -VF(x(y)) V 2 h(x(y)) + £ ti,V 2 f,(x(y)) (VF(x(y)))' 


s= 1 


5 Appendix: Code 

5.1 dual.R 

library (MASS) 
library (lsei) 
library (numDeriv) 

source ( "nnnewton.R" ) 
source (" qpqc.R" ) 


5.2 nnnewton.R 

nnnewton <- 

function (yold, 

fnewton, 
itmax = 100, 
eps = le-10, 
verbose = FALSE, 

...) { 
itel <- 1 
fold <- -Inf 

hfunc <- function (step, yold, ybar, ...) { 
z <- yold + step * (ybar - yold) 
fval <- fnewton (z, ...) 
return (fval$h) 

> 

repeat { 

fval <- fnewton (yold, ...) 
fnew <- fval$h 
ybar <- 

pnnqp (eps * diag(length (yold)) -fval$d2h, -drop(fval$dh - fval$d2h 7*7 yold))$: 
sval <- fnewton (ybar, ...) 
sd <- sum ((ybar - yold) * sval$dh) 
if (sd >= 0) step <- 1 
else { 
step <- 

optimize ( 
hfunc, 
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C (0, 1), 

yold = yold, 
ybar = ybar, 
maximum = TRUE, 

)$maximum 

} 

ynew <- yold + step * (ybar - yold) 
if (verbose) 
cat ( 

"Iteration: ", 

formate (itel, width = 3, format = "d"), 
"fold: ", 
formate ( 
fold, 

digits = 8, 
width = 12, 
format = "f" 

), 

"fnew: ", 
formate ( 
fval$h, 
digits = 8, 
width = 12, 
format = "f" 

), 

"step: ", 
formate ( 
step, 

digits = 8, 
width = 12, 
format = "f" 

), 

"\n" 

) 

if ((itel == itmax) II ((fval$h - fold) < eps)) 
break 

itel <- itel + 1 
yold <- ynew 
fold <- fnew 

> 

return (list ( 
solution = yold, 
value = fval$h. 
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gradient = fval$dh, 
hessian = fval$d2h, 
itel = itel 


} 


)) 


5.3 qpqc.R 

qpqc <- 

function (yold, 

a, 

b, 

c, 

itmax = 100, 
eps = le-10, 
analytic = TRUE, 
verbose = FALSE) { 

hfunc <- ifelse (analytic, hfgh, hfghnum) 
hval <- 

nnnewton ( 
yold, 
hfunc, 
a = a, 
b = b, 
c = c, 
eps = eps, 
verbose = verbose, 
itmax = itmax 

) 

fval <- primal (hval$solution, a, b, c) 
return ( 
list ( 

h = hval$value, 
f = fval$value, 
xmin = fval$solution, 
multipliers = hval$solution, 
constraints = hval$gradient 

) 

) 

} 

qfunc <- function (x, a, b, c) { 

return (c + sum (b * x) + sum (x * (a °/ 0 *°/o x)) / 2) 

> 


hfgh <- function (x, a, b, c) { 
m <- length (x) 
n <- nrow (a[, , 1] ) 
asum <- a[, , 1] 
bsum <- b [, 1] 
csum <- c[1] 
for (j in l:m) { 

asum <- asum + x[j] * a[, , j + 1] 

bsum <- bsum + x[j] * b[, j + 1] 

csum <- csum + x[j] * c [j + 1] 

} 


vinv <- solve (asum) 

xmin <- -drop (vinv bsum) 

h <- csum + sum (bsum * xmin) / 2 

dh <- rep (0, m) 

dg <- matrix (0, n, m) 

for (j in l:m) { 

dh [j ] <- qfunc (xmin, a[, , j + 1] , b[, j + 1], c [ j + 1] ) 
dg[, j] <- drop (b[, j + 1] + a[, , j + 1] xmin) 


d2h <- -crossprod (dg, vinv dg) 
return (list ( 
h = h, 


dh = dh, 


d2h = d2h 


} 


)) 


hfghnum <- function (x, a, b, c) { 
hfunc <- function (x) { 

return (hfgh (x, a, b, c)$h) 

} 

return (list (h = hfunc (x), dh = grad (hfunc, x), d2h = hessian (hfunc, x))) 


primal <- function (x, a, b, c) { 
m <- length (x) 
asum <- a[, , 1] 
bsum <- b [, 1] 
csum <- c[1] 
for (j in l:m) { 

asum <- asum + x[j] *a[, , j + 1] 

bsum <- bsum + x[j] *b[, j + 1] 

csum <- csum + x[j] * c [j + 1] 
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} 

xmin <- -drop (solve (asum, bsum)) 

fmin <- qfunc (xmin, a[, , 1], b[, 1], c[l]) 

return (list (solution = xmin, value = fmin)) 

> 
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